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Precision  Orbit  Determination  (OD)  is  often  critical  for  successful  satellite  operations 
supporting  a  wide  variety  of  missions.  Directional  or  angles  only  measurements  of  a  satel¬ 
lite  are  typically  represented  in  spherical  coordinates  on  the  observer’s  celestial  sphere 
(e.g.  azimuth  and  elevation  or  right  ascension  and  declination).  Computing  residuals  in 
these  angular  coordinates  during  the  OD  process  can  introduce  errors  in  the  same  way  an 
equirectangular  projection  distorts  both  distance  and  direction  on  a  map.  One  technique 
is  to  weight  the  azimuth  residuals  by  the  cosine  of  the  elevation,  as  in  a  sinusoidal  projec¬ 
tion.  While  this  technique  preserves  the  length  of  every  parallel  it  still  induces  distortion 
in  direction  and  distance.  Therefore,  it  is  proposed  to  use  the  angular  distance  between 
the  computed  and  observed  locations  as  the  residual.  This  is  similar  to  using  an  azimuthal 
equidistant  projection  with  the  observed  location  at  the  center.  It  is  shown  that  this  tech¬ 
nique  removes  distortion  present  in  the  other  two  representations.  The  three  techniques 
are  then  compared  experimentally  for  a  geostationary  and  a  low  Earth  orbit  satellite  using 
simulated  data  to  evaluate  their  differences.  It  is  shown  that  using  angular  distance  as 
the  residual  decreases  the  number  iterations  required  for  convergence  and  allows  the  OD 
process  to  more  closely  fit  the  observed  data  when  there  are  observations  near  the  pole  of 
the  spherical  coordinate  representation. 

I.  Introduction 

Accurate  Orbit  Determination  (OD)  is  often  critical  for  successful  satellite  operations  supporting  a  wide 
variety  of  missions.  Precision  OD  involves  accurately  modeling  all  satellite  observations,  including  angular 
or  directional  observations  such  as  those  obtained  from  astrometry.  Standard  methods  for  OD,  the  Kalman 
Filter  (KF)  and  Weighted  Least  Squares  (WLS),  assume  a  Gaussian  distribution  for  the  observations  and  then 
attempt  to  find  the  Maximum  Likelihood  Estimate  (MLE).  Therefore  the  choice  of  observation  representation 
has  a  significant  effect  on  the  size  and  shape  of  the  likelihood  function  which  it  turn  affects  the  MLE.  Other 
work  on  representation  in  OD  has  focused  on  the  estimated  state  where  Junkins  et.  al.1  found  the  choice  of 
state  representation  has  a  significant  effect  on  the  shape  and  characteristics  of  the  estimated  covariance. 

Standard  OD  references  suggest  representing  directional  measurement  as  azimuth  and  elevation;  right 
ascension  and  declination;  North-South  and  East- West  angles;  or  some  other  representation  of  direction  in 
spherical  coordinates. 2-6  An  orbit  determination  manual  by  Kuga  and  Gill'  describes  a  process  of  converting 
angular  measurements  to  a  unit  vector  to  apply  refraction  and  then  converting  back  to  spherical  coordinates 
for  computing  the  residuals.  Similarly  the  Goddard  Trajectory  Determination  System  (GTDS)  mathematical 
theory8  uses  spherical  coordinates  representation  for  direction  measurements.  Several  recently  published 
works  on  OD  compute  residuals  in  spherical  coordinates  as  well.9, 10  The  Naval  Research  Laboratory’s  (NRL) 
Orbit  Covariance  Estimation  and  Analysis  (OCEAN)  OD  application  weights  azimuth  residuals  by  the  cosine 
of  the  corresponding  elevation  angle  to  account  for  the  convergence  of  meridians  near  the  poles.11  None  of 
these  references  suggest  using  angular  distance  on  the  celestial  sphere  as  the  residual  as  is  recommended  in 
this  paper. 

Choosing  a  representation  for  directional  measurements  is  equivalent  to  choosing  a  transformation  from 
location  on  a  sphere  to  location  in  a  plane  which  has  a  long  history  in  map  making.  Using  a  point’s  spherical 
coordinates  as  its  planar  coordinates  is  known  as  the  equirectangular  or  plate  carree  projection.  While  simple 
to  construct,  it  is  problematic  because  of  its  distortion  in  angle  and  distance,  especially  near  the  poles.12 
Weighting  longitude  by  the  cosine  of  latitude  is  the  sinusoidal  projection  which  has  “no  distortion  along  the 
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Equator  and  central  meridian,  but  distortion  becomes  pronounced  near  the  outer  meridians,  especially  in 
the  polar  regions.” 12  The  azimuthal  equidistant  projection,  a  projection  that  is  often  used  for  air  navigation, 
preserves  distance  and  direction  from  the  central  point.12  The  distortion  caused  by  the  chosen  projection 
changes  how  the  OD  process  computes  distance  (residual)  and  direction  (gradient)  between  observed  and 
computed  measurements. 

In  the  following  section,  three  direction  observation  representations  are  defined  and  compared  analytically. 
The  three  representations  are  spherical  coordinates,  weighted  spherical  coordinates,  and  a  unit  vector  based 
approach.  Then  Section  III  describes  a  method  for  comparing  the  representations  using  simulated  data  for  a 
Low  Earth  Orbit  (LEO)  and  a  Geostationary  Earth  Orbit  (GEO)  satellite.  Section  IV  discusses  the  results 
from  the  numerical  simulations  and  finally  Section  V  presents  the  authors’  conclusions. 

II.  Direction  Representation 

The  choice  of  direction  representation  affects  not  only  the  size  of  the  residual  but  also  the  shape  of  the 
assumed  probability  distribution  since  in  WLS  OD  or  KF  OD  the  observations  are  assumed  to  be  samples 
of  a  Normal  distribution.  For  the  OD  solution  to  be  the  Maximum  Likelihood  Estimate  (MLE),  the  cost 
function  must  have  the  same  shape  as  the  probability  distribution  of  the  observations.  For  this  paper  it  is 
assumed  observations  are  Normally  distributed  about  the  true  point  on  the  celestial  sphere,  the  distribution 
is  symmetric  about  the  true  point  and  the  distribution  is  independent  of  location  on  the  celestial  sphere. 

All  angular  direction  representations  incur  some  error  because  the  domain  of  the  Normal  distribution  is 
infinite  but  the  range  of  any  angular  measurement  is  at  most  2n,  though  this  limitation  is  minimal  when 
the  limits  are  many  standard  deviations  away  from  the  evaluated  point.  A  potentially  larger  concern  is 
the  local  distortion  in  shape  and  distance  caused  by  the  choice  of  projection  from  the  sphere  to  the  plane. 
The  following  discussions  examine  three  different  direction  representations  that  correspond  to  the  the  plate 
carree,  sinusoidal  and  azimuthal  equidistant  projections. 

Spherical  Coordinates 

The  first  approach  is  directly  using  spherical  coordinates  where  direction  is  represented  as  a  pair  of  angles 
such  as  azimuth  and  elevation;  longitude  and  latitude;  or  right  ascension  and  declination.  The  following 
discussion  uses  azimuth  and  elevation  though  mathematically  these  systems  are  equivalent,  differing  only  in 
the  choice  of  reference  planes.  Residuals  formed  in  spherical  coordinates  take  the  form 


ra  =  ac-  a0  =  Act  (1) 

r£  =  ec  -  e0  =  Ae  (2) 

where  r  is  the  residual,  a  is  the  azimuth  angle,  £  is  the  elevation  angle,  subscript  o  denotes  observed  and 
subscript  c  denotes  computed.  The  corresponding  cost  function  to  be  minimized  in  the  OD  process  is 

Cs=rl+r2s  =  (  Aa)2  +  (Ae)2.  (3) 

This  cost  function  is  analogous  to  using  an  equirectangular  projection.  Though  effective  near  zero  elevation, 
the  azimuth  residuals  can  be  large  for  small  differences  in  direction  near  the  poles  (e  =  90°). 

Weighted  Spherical  Coordinates 

The  second  approach  is  to  preserve  area  in  mapping  the  celestial  sphere  to  the  plane  by  weighting  the  azimuth 
residual  by  the  length  of  the  corresponding  small  circle  of  constant  elevation, 

Cw  =  (Act  cos  eG)2  +  (Ae)2  (4) 

which  is  similar  to  the  sinusoidal  projection.  While  this  projection  preserves  area,  it  does  not,  in  general, 
preserve  distance  or  direction  between  points.  This  has  the  effect  of  de-weighting  azimuth  observations  near 
the  poles.  In  the  extreme  case  where  there  are  many  observations  at  high  elevation  this  representation  can 
reduce  the  observation  set  by  half  by  eliminating  the  azimuth  observations. 
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Unit  Vectors 

A  third,  alternative  representation  can  be  derived  from  using  unit  vectors  to  represent  directions.  The  angle 
between  the  observed  and  computed  directions  can  be  used  as  the  residual. 

re  =  9  (5) 

The  angle  between  two  vectors  can  be  computed  either  from  the  dot  product 

cos0  =  uc  •  u0  (6) 

or  the  cross  product 

sin  0  =  | uc  x  uQ  |  (7) 

where  u  is  a  unit  vector.  Though  mathematically  equivalent,  Equation  7  is  preferred  for  computer  imple¬ 
mentation  because  the  dot  product  is  insensitive  to  changes  when  6  is  small.  The  unit  vector  derived  cost 
function  is 

cu  =  e2.  (8) 

While  Equation  8  effectively  measures  the  angular  distance  between  observed  and  computed  points  it 
only  provides  one  equation  per  direction  observation.  A  second  angular  equation  can  be  added  to  keep  two 
independent  equations  per  direction  observation  by  measuring  the  angle  out  of  the  uc,  uQ  plane.  First  define 
the  normal  to  the  plane 

u„  =  uc  x  u0  (9) 

then 

sin  (j>  —  uc  ■  u„  (10) 

defines  the  out  of  plane  angle  and  =  (f>  is  the  residual.  When  differentiating  Equation  10  u„  is  treated  as 
a  constant  since  it  is  defining  a  new  frame.  The  residual  does  not  contribute  to  the  cost  or  its  derivative 
because  it  is  always  zero.  In  this  way  the  unit  vector  approach  can  be  viewed  as  defining  a  new  system  of 
spherical  coordinates  for  every  computed  point  with  uc  and  uQ  defining  the  equatorial  plane. 

Comparison 

In  general,  the  three  approaches  produce  different  results.  However,  under  certain  assumptions  Cu  is  shown 
to  be  equivalent  to  Cw  and  Cw  is  shown  to  be  equivalent  to  Cs.  This  result  is  used  Section  III  to  design 
experiments  that  show  the  differences  and  similarities  between  the  methods. 

For  small  values  of  9 ,  Aa  and  Ae  Cu  is  equivalent  to  Cw.  Using  the  small  angle  approximation  and 
Equation  7 

9  —  |  uc  x  uG  | .  (11) 

Expressing  each  unit  vector  in  spherical  coordinates  and  taking  the  magnitude  of  the  cross  product  gives 

Cu  =  92  =  (sacs£0c£c  —  sa0s£cc£0)2  +  (sec  ca0  c£a  —  s£0cacc£0)2  +  s2  (ac  —  a0)  c2  ecc2  eG  (12) 

where  s  and  c  denote  sine  and  cosine  respectively.  Substituting  ac  =  a0  +  A  a,  £c  =  £a  +  As,  and  ignoring 
term  of  order  greater  than  two  results  in 

Cu  =  (Aa)2  cos2  £0  +  (Ae)2  (13) 

which  is  Equation  4.  Therefore,  it  is  also  true  that  the  derivatives  Cu  are  equivalent  to  the  derivatives  of  Cw 
when  0,  Aa,  and  Ae  are  all  small. 

For  small  values  of  eQ  Cw  is  equivalent  to  Cs .  Using  the  approximation  that  cos  eD  «  1 

Cw  =  (Aa)2  +  (Ae)2 

which  is  the  definition  of  Cs  in  Equation  3.  Therefore,  when  9 ,  Aa,  Ae,  and  e0  are  all  small  the  three 
methods  are  equivalent. 

While  spherical  and  unit  vector  based  representations  are  the  same  for  small  angles  they  produce  different 
results  when  any  one  of  three  angles  ( 9 ,  Aa,  or  Ae)  is  large.  For  large  angles  Cu  still  accurately  measures 
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distance  on  the  celestial  sphere  because  that  is  its  definition,  but  Cw  and  Cs  do  not  because  they  are  not 
equivalent  to  Cu  for  large  angles.  Large  residuals  may  occur  when  the  initial  guess  of  the  satellite’s  orbit  is 
far  from  its  actual  orbit,  as  is  common  in  initial  OD  scenarios.  The  distortion  caused  by  using  a  spherical 
representation  is  hypothesized  to  cause  slow  convergence  or  cause  divergence. 

When  the  initial  guess  is  accurate,  9  and  Ae  are  small  but  Aa  may  still  be  large  when  there  are  observa¬ 
tions  near  the  pole  of  the  spherical  representation.  An  example  is  illustrated  in  Figure  1  where  9  <  5°  only 
implies  Aa  <  90°  which  is  quite  large.  From  the  figure  it  is  apparent  that  there  are  significant  differences 
between  the  cost  functions,  i.e.  assumed  error  distributions,  imposed  by  the  different  representations.  Cw 
favors  lower  elevations  than  Cu  and  Cs  is  much  thinner  in  azimuth,  hence  emphasizing  azimuth  measure¬ 
ments  more  than  the  other  two  cost  functions.  Since  the  gradient  is  normal  to  the  level  curves  shown  in 
Figure  1  it  is  also  evident  there  are  significant  differences  in  the  derivatives.  For  example,  near  the  cusp 
of  Cw  at  e  =  90°  the  gradient  is  normal  to  the  direct  great  circle  path  to  uc.  As  well  as  causing  slow 
convergence,  the  differences  in  the  cost  function  may  cause  the  OD  process  to  converge  to  different  solutions. 
It  is  hypothesized  that  using  a  unit  vector  representation  leads  to  more  accurate  OD  than  using  a  spherical 
coordinate  representation  when  there  are  observations  very  close  to  the  pole. 

Cost  Function  Distortion  Near  The  Pole 


+85°  +85° 


Figure  1:  Level  curves  of  the  three  cost  functions  near  the  poles  showing  an  example  of  the  distortion  caused 
by  using  different  direction  representations.  The  figure  uses  an  azimuthal  equidistant  projection  centered  on 
the  observed  direction.  Each  level  curve  is  plotted  for  the  value  C  =  (5°)2. 


III.  Method 

The  three  direction  representations  are  tested  with  simulated  data  using  a  NRL  WLS  OD  software 
program  called  State  Estimation  Application,  or  SEA.  This  tool  is  written  in  Java  using  the  Orekita  space 
dynamics  library  and  Hipparchusb  mathematics  library.  SEA  is  used  at  NRL  for  OD  research  and  is  compared 
to  NRL’s  operational  OD  application,  OCEAN,  in  a  previous  paper.13 

Each  representation  is  tested  for  a  Low  Earth  Orbit  (LEO)  and  a  Geostationary  Earth  Orbit  (GEO) 
satellite  with  the  initial  orbital  parameters  shown  in  Table  1.  Different  ground  sites  are  used  for  the  LEO  and 


ahttps:/ /orekit.org/ 
bhttps:/ / www. hipparchus.org/ 
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Table  1:  Orbital  Parameters  for  the  LEO  and  GEO  satellite  Table  2.  Ground  Stations  used  for  LE0  and 
_ _ _  GEO  satellites 


Parameter 

LEO 

GEO 

Semimajor  Axis 
Eccentricity 
Inclination 
Argument  of  Perigee 
Lon.  of  Ascending  Node 
True  Anomaly 

6878.14  km 

0 

45  deg 

0  deg 

235.929  deg 

0  deg 

42,166.3  km 

0 

0.01  deg 

0  deg 
30.018  deg 
360  deg 

GS-LEO 

GS-GEO 

Latitude 

Longitude 

Altitude 

38.9  N  deg 
77.0  W  deg 

0  km 

0.0  deg 
30.0  deg 

0  km 

GEO  test  cases  as  shown  in  Table  2.  The  LEO  test  cases  are  designed  to  evaluate  the  different  representations 
when  all  three  are  expected  to  have  similar  performance  since  for  LEO  satellite  passes  most  of  the  data  has 
low  elevation  angles.  The  GEO  satellite’s  orbit  is  chosen  to  evaluate  an  extreme  case  where  all  observations 
are  close  to  the  pole. 

Truth  ephemerides  are  produced  in  Systems  Took  Kit  (STK)  using  its  High-Precision  Orbit  Propagator 
(HPOP).  STK  is  used  instead  of  SEA  to  produce  the  truth  ephemeris  so  that  different  implementation  of 
the  force  force  and  measurement  models  are  used  during  OD.  Both  truth  ephemerides  include  the  WGS84 
-  EGM96  gravity  model,  ocean  tides,  and  Sun  &  Moon  third  body  gravity.  The  GEO  ephemeris  includes 
solar  radiation  pressure.  The  same  force  models  were  included  in  SEA  during  OD  to  maintain  consistency. 


(a) 


Error  Distribution  Near  The  Pole 

+89.6°  +89.8°  -135°  180°  +135°  +89.8°  +89.6° 


Figure  2:  Observation  error  distribution  with  samples  at  two  example  points.  Each  plotted  circle  has  a 
radius  of  1,  2,  and  3  times  the  standard  deviation.  The  plot  on  the  left  is  centered  at  0°  azimuth  0°  elevation 
and  the  plot  on  the  right  is  centered  at  0°  azimuth  89.7°  elevation. 


The  truth  ephemeris  is  used  to  create  simulated  input  observations  for  SEA.  Noise  is  added  to  the  direction 
observations  such  that  the  distribution  is  symmetric  about  the  observed  direction  and  is  independent  of  the 
observed  direction’s  location  on  the  celestial  sphere.  Noise  is  Normally  distributed  with  a  standard  deviation 
of  0.1  degree.  The  distribution  is  shown  for  two  example  points  in  Figure  2. 

The  LEO  test  cases  have  observations  that  consist  of  several  passes  over  a  given  day  as  shown  in  Table  3. 
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Table  3:  LEO  Observations 


Pass  Number 

Start  Time  (UTC) 

Duration  (min) 

Max  Elevation  (deg) 

Used  in  LEO  Test  Case  2 

P-1 

00:11 

9.68 

79.43 

V 

P-2 

01:50 

9.29 

36.51 

P-3 

03:30 

9.22 

35.05 

P-4 

05:09 

9.61 

88.38 

V 

P-5 

06:48 

7.57 

16.57 

P-6 

22:09 

7.33 

14.97 

P-7 

23:46 

9.68 

80.63 

V 

All  seven  of  the  passes  are  used  in  the  first  LEO  test  case,  while  only  the  passes  that  have  a  high  maximum 
elevation  angle  were  are  in  the  second  LEO  test  case.  The  LEO  ground  station  is  set  to  have  an  elevation 
mask  of  5  degrees.  The  observations  are  generated  every  10  seconds. 

The  GEO  test  cases  have  observations  spanning  an  entire  day  with  a  time  step  of  30  seconds.  The  first 
GEO  test  case  uses  every  observation  for  a  total  of  2880  data  points.  The  second  GEO  test  case  uses  the  same 
set  of  observations,  but  has  several  black-out  periods  of  1-2  hours  to  evaluate  the  different  representations 
performance  with  a  reduced  data  set.  The  total  number  of  observations  used  in  the  this  case  is  1500. 

Once  generated,  the  observations  are  processed  by  SEA  using  a  separate  measurement  model  for  each 
direction  representation.  In  each  case  the  weights  are  set  appropriately  for  a  standard  deviation  of  0.1  degrees. 
SEA  uses  a  Levenburg-Marquardt  optimizer  and  the  convergence  criteria  is  a  change  in  the  weighted  RMS 
residuals  of  10-6  or  less. 

Table  4:  The  magnitude  of  the  initial  condition  offset  from  truth. 


Case 

IC-1 

IC-2 

IC-3 

IC-4 

LEO 

GEO 

214.0  km 

153.9  km 

420.1  km 

35.9  km 

886.8  km 

10.2  km 

395.4  km 

8.1  km 

GEO  Initial  Conditions 
N 


(a)  Snapshot  of  LEO  ICs.  The  truth  orbit  is  located  on  the  (b)  View  of  the  GEO  ICs  as  seen  from  the  ground  station, 
far  right. 

Figure  3:  Initial  Conditions  for  OD. 


A  single  run  evaluates  a  single  direction  representation  with  a  set  of  observations  and  a  given  initial 
condition.  For  each  initial  condition  three  runs  are  performed,  one  with  each  direction  representation.  A 
test  case  consists  of  using  four  different  initial  conditions  with  the  same  observation  set.  The  initial  conditions 
are  shown  in  Table  4  and  graphically  in  Figure  3.  The  initial  conditions  are  chosen  to  evaluate  the  different 
representations  for  different  magnitudes  of  the  initial  resdiuals.  For  each  satellite  one  test  case  is  run  with 
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all  the  observations  and  then  a  second  test  case  is  run  with  the  reduced  set  of  observations,  for  a  total  of 
four  test  cases. 


IV.  Results 

The  direction  representations  are  evaluated  based  on  two  criteria:  the  difference  between  estimated  and 
truth  ephemerides  during  the  fit  span,  and  the  number  of  iterations  to  convergence.  The  ephemeris  difference 
is  computed  as  the  RMS  of  the  magnitude  of  the  position  difference  at  every  point  during  the  fit  span. 

LEO  Test  Cases 


LEO  Test  Case  1 


1400 


1200 


12  3  4 


Axis  Title 


■  Unit  Vector 

■  Spherical  Coordinates 

Weighted  Spherical 
Coordinates 


Figure  4:  LEO  test  case  1  ephemeris  difference. 


LEO  Test  Case  2 


Initial  Condition 


■  Unit  Vector 

■  Spherical  Coordinates 

Weighted  Spherical 
Coordinates 


Figure  5:  LEO  test  case  2  ephemeris  difference. 


Ephemeris  differences  for  the  LEO  test  cases  are  shown  in  Figures  4  and  5.  In  all  cases  the  weighted 
spherical  coordinates  representation  and  the  unit  vector  representation  converge  to  nearly  the  same  solution, 
differing  by  up  to  0.2%  in  the  estimated  ephemeris’s  distance  from  the  truth  ephemeris.  The  spherical  coordi¬ 
nates  representation  always  converges  to  a  solution  farther  from  the  truth  than  the  other  two  representations. 
The  spherical  coordinates  representation  has  a  RMS  at  least  1.8  times  greater  than  the  other  representations 
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Figure  6:  The  number  of  iterations  until  convergence  for  the  first  LEO  test  case. 
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Figure  7:  The  number  of  iterations  until  convergence  for  the  second  LEO  test  case. 


in  LEO  test  case  1  and  at  least  1.3  times  greater  in  LEO  test  case  2.  When  starting  from  IC-3,  the  furthest 
initial  condition,  each  representation  converges  to  a  different  solution  than  when  started  from  the  other  three 
initial  conditions.  This  suggests  multiple  local  optima  exist  and  that  the  initial  condition  determines  which 
is  selected. 

The  number  of  iterations  for  each  method  are  shown  in  Figures  6  and  7.  Spherical  coordinates  repre¬ 
sentation  uses  more  iteration  than  the  other  two  representations  in  every  case.  Unit  vector  representation 
converges  in  fewer  iterations  than  weighted  spherical  coordinates  representation  in  six  out  of  the  eight  cases 
examined.  Reducing  the  amount  of  data  does  not  seem  to  have  a  significant  effect  on  the  number  of  iterations 
until  convergence.  Another  trend  is  that  all  methods  converge  the  quickest  using  IC-1,  which  is  the  closest 
to  the  truth,  and  slowest  using  IC-3,  which  is  the  farthest. 

GEO  Test  Cases 

Compared  to  the  other  methods  the  unit  vector  representation  produces  an  estimate  with  the  smallest 
ephemeris  difference  as  shown  in  Figures  8  and  9.  The  RMS  ephemeris  difference  for  the  unit  vector  repre¬ 
sentation  is  3.418  km  in  GEO  test  case  1  and  3.187  km  in  GEO  test  case  2  for  all  initial  conditions  tested. 
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Figure  8:  GEO  test  case  1  epliemeris  difference. 
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Figure  9:  GEO  test  case  2  epliemeris  difference. 
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Figure  10:  The  truth  orbit  and  three  orbit  solutions  each  determined  using 
a  different  direction  representation  as  seen  from  the  GS  ground  station  point 
of  view. 
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Figure  11:  The  number  of  iterations  until  convergence  for  the  first  GEO  test  case. 
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Figure  12:  The  number  of  iterations  until  convergence  for  the  second  GEO  test  case. 


This  implies  the  unit  vector  representation  converges  to  the  same  solution  regardless  of  initial  condition.  The 
unit  vector  representation  achieves  a  lower  RMS  ephemeris  difference  than  the  other  two  representations  by  a 
factor  of  at  least  2.6  in  all  test  cases.  The  greatest  difference  appears  in  test  case  1  IC-2  where  the  ephemeris 
difference  for  the  unit  vector  representation  is  47.3  times  smaller  than  the  spherical  coordinates  representa¬ 
tion.  Additionally,  spherical  coordinates  representation  leads  to  significant  variation  in  the  solution  while 
the  other  two  representations  consistently  converge  to  the  same  solution.  Figure  10  shows  the  truth  orbit 
of  the  GEO  satellite  as  well  as  solutions  using  the  three  different  methods  from  the  ground  station’s  point 
of  view,  illustrating  the  differences  in  solution  accuracy.  In  the  figure  the  estimate  from  the  unit  vector 
representation  is  almost  indistinguishable  from  the  truth  ephemeris.  The  solutions  shown  are  from  GEO 
test  case  1,  IC-3. 

The  unit  vector  method  converges  within  3  iterations  for  each  IC  in  both  GEO  test  cases  which  is  fewer 
than  then  other  methods  as  shown  in  Figures  11  and  12.  Even  with  fewer  observations,  the  unit  vector 
method  was  able  to  converge  to  a  solution  using  the  same  number  of  iterations.  The  two  other  methods  took 
longer  to  converge  when  there  were  fewer  observations  within  the  fit  span.  In  all  except  for  GEO  test  case 
2  IC-1,  the  weighted  spherical  coordinates  representation  took  a  few  iterations  longer  to  converge  than  the 
unweighted  spherical  coordinates  representation. 
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The  first  hypothesis,  that  the  unit  vector  representation  converges  quicker  when  the  errors  in  the  initial 
condition  are  large,  is  confirmed.  In  the  GEO  test  cases  IC-1  is  the  farthest  initial  condition  and  the  unit 
vector  representation  converges  quicker  than  the  other  methods  by  at  least  13  iterations.  For  the  LEO  test 
cases  IC-3  is  the  farthest  and  the  unit  vector  representation  converges  quicker  than  the  other  methods  by 
at  least  three  iterations.  When  the  errors  in  the  initial  condition  are  small,  as  for  LEO  IC-1,  unit  vector 
representation  uses  the  same  number  of  iterations  or  one  more  iteration  than  weighted  spherical  coordinates 
representation. 

The  second  hypothesis,  that  using  the  unit  vector  representation  produces  a  more  accurate  solution 
in  the  presence  of  high  elevation  observations,  is  confirmed  in  the  GEO  case  where  all  observations  are 
above  89  degrees  elevation.  The  unit  vector  approach  improves  solution  accuracy  by  more  than  an  order  of 
magnitude  while  reducing  the  number  of  iteration  by  more  than  a  factor  of  five.  This  result  is  partly  due 
to  the  distortion  caused  by  weighted  and  un-weighted  spherical  coordinate  representation  at  high  elevation 
angles.  This  result  is  also  partly  due  to  the  assumed  probability  distribution  in  elevation  extending  beyond 
the  physical  bounds  of  the  elevation  angle  when  using  a  weighted  or  un-weighted  spherical  coordinates 
representation.  While  a  physics  based  elevation  distribution  has  zero  probability  of  £  >  90°  the  assumed 
distribution  places  significant  probability  on  non-physical  elevation  angles,  leading  to  a  less  accurate  solution. 
In  the  weighted  spherical  coordinates  representation  the  azimuth  observations  are  effectively  deleted  from  the 
dataset  by  assigning  them  a  near  zero  weight  which  reduces  the  dataset  by  half.  In  the  unweighted  spherical 
coordinates  representation  far  too  much  emphasis  is  placed  on  azimuth  observations  which  appears  to  cause 
the  OD  process  to  fit  some  of  the  noise  in  the  observations.  Both  spherical  coordinates  representations  are 
unable  to  average  the  remaining  elevation  observations  to  90°  because  the  average  elevation  angle  in  the 
input  data  is  89.9°. 

The  LEO  test  cases  include  observations  for  a  range  of  elevation  angles  between  5  and  88  degrees.  For 
low  elevation  angles,  a  majority  of  the  data,  a  small  6  (angle  between  observed  and  computed  points)  implies 
small  Aa  and  As  so  it  is  not  surprising  that  the  unit  vector  and  weighted  azimuth  elevation  representations 
converge  to  the  same  solution  as  they  have  the  same  cost  function  in  these  cases.  While  this  finding  provides 
no  direct  evidence  for  or  against  the  second  hypothesis  it  does  indicate  indicate  the  limits  of  where  it  applies. 
Namely  that  a  larger  percentage  of  higher  elevation  observations  is  needed  for  unit  vector  representation 
to  produce  a  more  accurate  solution  than  weighted  spherical  coordinates  representation.  The  unweighted 
spherical  coordinates  representation  produced  a  less  accurate  ephemeris  than  the  other  two  representations 
in  every  case.  For  very  low  elevation  observations  (cose  «  1)  all  three  representations  are  equivalent  but  it 
seems  that  by  overweighting  azimuth  observations  at  higher  elevations  the  unweighted  spherical  coordinates 
representation  uses  a  significantly  different  cost  function  which  leads  to  a  less  accurate  solution. 

V.  Conclusion 

In  conclusion,  the  unit  vector  representation  outperforms  weighted  and  unweighted  spherical  coordinates 
representation  when  there  are  a  large  number  of  observations  very  close  to  the  pole.  It  is  both  more  accurate 
in  the  estimated  solution  and  able  to  converge  to  that  solution  in  fewer  iterations.  When  there  are  few 
observations  near  the  pole,  unit  vector  representation  and  weighted  spherical  coordinates  representation  are 
equivalent  in  the  accuracy  achieved  while  using  unweighted  spherical  coordinates  still  leads  to  less  accurate 
solutions.  The  improvement  due  to  the  unit  vector  representation  is  due  to  more  accurately  fitting  the  error 
distribution  of  the  observations.  Therefore,  the  authors  recommend  using  the  unit  vector  representation 
proposed  in  this  paper  when  the  error  in  the  observations  is  Normally  and  symmetrically  distributed. 

These  results  are  obtained  from  a  limited  set  of  experimental  runs  with  simulated  data;  examining  a 
wider  variety  of  cases  would  better  confirm  the  results  as  well  find  the  tipping  point  where  the  unit  vector 
representation  outperforms  the  other  representations.  Further  work  could  be  performed  to  characterize  the 
effect  of  direction  measurement  representation  on  the  estimated  covariance.  Additional  future  work  could 
tailor  a  direction  measurement  representation  to  an  empirical  distribution  derived  from  actual  observations. 

Finally,  by  providing  quicker  and  more  accurate  OD  solutions  the  new  unit  vector  based  direction  repre¬ 
sentation  can  contribute  to  the  success  of  a  wide  variety  of  satellite  missions. 
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